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^ ; ABSTRACT 
Oh! 

. We consider the patterns formed by small rod-like objects advected by a random flow 

i in two dimensions. An exact solution indicates that their direction field is non-singular. 

However, we find from simulations that the direction field of the rods does appear to 
exhibit singularities. First, 'scar lines' emerge where the rods abruptly change direction 
' by TT. Later, these scar lines become so narrow that they 'heal over' and disappear, but 

^ ■ their ends remain as point singularities, which are of the same type as those seen in 

fingerprints. We give a theoretical explanation for these observations. 



^ . 1. Introduction 

o 

' We consider the motion of small rod-like particles suspended in a moving fiuid. The suspended 

ff^ ■ particles align with their neighbours in a manner determined by the strain-rate of the fiow. In a 

Q> ■ turbulent or randomly moving fiuid the direction vector field of the rods forms complex textures, 

^ . illustrated by figure [H We concentrate on two-dimensional textures, because it is hard to observe 

, the direction field in three dimensions. Also, we confine attention to the case of incompressible 

flow, which is easiest to analyse and which is easily realised experimentally (by using a suspension 
of rod-like particles in a fllm of water floating upon a denser fluid which is randomly stirred). The 



>< 



?H ' results are of quite general interest, because any asymmetric particles will have a preferred direction 

determined by the history of the strain tensor of the field along the trajectory of the particle. 

Suspen sions of small anisotropic particles cafled rheoscopic fluids are often used for flow vi- 



sualisation (jMatisse fc Gormanlll984l : ISava§lll985l : iGauthier. Gondoret Rabaudlll998l ). This uses 



the principle that the intensity of scattering of light from a localised source will depend upon the 
orientation of the suspended particles. T he information in this visu alisation can be enhanced by 



using light sources with different colours (iThoroddson &: Bauerjll999l ). In this paper we show how 



the colours might be used to reveal information about the topology of the textures formed by the 
rheoscopic fluid. In flgure [2] we demonstrate the potential of this approach for rod-like particles. 
For illustrative purposes, we assume that the intensity of scattering from a rod at angle 9 from a 
source at angle (p (relative to a line perpendicular to the rod) is proportional to cos^(^ — cp) (this 
approximation can be justifled when the rods are short compared to the wavelength of the light). 
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Fig. 1. — Simulations of the orientations of rods advected by a random flow in two dimensions. 
These appear to show singularities which are analogous to those occurring in fingerprint patterns, 
illustrated in figure El The details of the simulations are specified in the appendix. 



Accordingly, in figure [2] we re-display the textures in figure [T] by plotting a colour C which is an 
admixture of the primary colours red, green and blue, denoted (i?, G, B): 

C = Rcos^{e) + Gcos^{e - 2^/3) + Bcos^{e - 47r/3) . (1) 

(In figure [2] the angle 9 of the rods is measured relative to the horizontal, with increasing in 
the anti-clockwise direction.) The physics of scattering or reflection from the rod- like particles is 
complex, but this illustration is indicative of what can be seen with different coloured light sources. 

The rod direction field is a non-oriented vector field in a two-dimensional space (by non- 
oriented, we mean that rod directions differing by vr are equivalent). In such a field we might 
expect to see point singularities of the direction fiel d of the type illustrated in figure [3l which 



are also present in fingerprint patterns (iHenrvl 1 19001 ) (where the patterns formed by ridges are 



another example of a non-oriented vector field in two dimensions). The actual textures that we 
observe in simulations do indeed have structures which resemble the core and delta singularities of 
fingerprints, as illustrated by the examples in figured! We shall argue that the principles underlying 
the structures visible in this picture are quite subtle, and that it is in fact surprising to see such 
singularities. We remark that the singularities are characterised by a topological invariant, termed 
the Poincare index, which is illustrated in figure HI and topological arguments will be central to the 
discussion. Singularities with a non-zero Poincare index could be detected using the visualisation 
technique illustrated in figure [2l by examining the colours along a closed path. If the colours cycle 
through all three primaries as the path is traversed, this path must contain a singularity of the rods 
directions. The sign if the Poincare index is determined by the order in which the primary colours 
cycle (R ^ G ^ B or R ^ B ^ G). 
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Fig. 2. — The rod textures shown in figure [T] colour-coded using equation ([T]), to illustrate how the 
textures can be visualised using coloured light sources. 



In section[2]below we give a simple derivation of the equation of motion for th e rods, and pre sent 
its general solution. Our equation of motion is a limiting case of that given by iJefFeryl (119221 ) for 
the motion of an ellipsoid of revolution in a viscous fluid at low Reynolds number. We remark that 
a solution of this equ ation has bee r i obtained fo r a simple shear flow, in which the ellipsoid exhibits 
a tumbling motion (IJefFeryl Il922l : ISavaa Il985l ) , and that several authors have discussed chaotic 
aspects of this tumbling motion i n more complex flow s (IShin &: Maxeyi Il99ll : iMallier &: Maxey 
199ll : ISzeri. Wiggins k. Leall Il99ll : IShin &: MaxevI 119971). Suspens ions of anisotropic particles can 
be used for visualisation of fluid flows (IMatisse &: GormanI Il984l'). and the interpre tation of the 
images produced by these fluids is discussed bv lGauthier. Gondoret &: RabaudI (jl998l ). In section [2] 
we present our solution of the equation of motion, in terms of the monodromy matrix of the flow. In 
this paper we present a solution which is specific to the limiting case of rod-like particles; the case 
of more general shapes and three-dimensional flows will be discussed in a later work. We also show 
that the vector field is asymptotic (at large times) to the vector field formed by the eigenvectors of 
the monodromy matrix corresponding to the larger eigenvalue. 

Section [3] discusses the extent to which the solution we obtain in section [2] can exhibit sin- 
gularities. We start by presenting an argument showing that the direction field cannot have any 
singularities. This implies that the Poincare index for any curve is zero, and is hard to reconcile 
with the appearance of figures (H [2j Throughout most of the plane the direction field of the rods is 
asymptotic to the eigenvector field of the monodromy matrix. However, the eigenvector field can 
have a non-zero Poincare index implying that the asymptotic correspondence between these vector 
fields breaks down somewhere. We show that it fails along certain lines, which we term scar lines, 
where the direction vector of the rods abruptly changes by vr. The scar line emerges and sharpens 
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Fig. 3. — The t extures illustr ated in figure [T] have similarities with fingerprints patterns, such as 
(a) (taken from iHenryl ([l900|))- Such patterns contain two elementary point singularities of non- 
oriented vector fields in two-dimensions: in fingerprint patterns these are known as the core (b) and 
the delta (c), marked by red and green circles respectively in (a). Examples of these singularities 
as they appear in rod textures are shown in (d) and (e) respectively. 



as the two vector fields asymptotically approach each other. As the scar line (illustrated in figure 
[5]) sharpens, there will be fewer rods which lie in the region where the direction diff'ers from the 
asymptotic value. The result is that the scar line disappears (see figure [6|). At the ends of the scar 
line there remains a point singularity of the type illustrated in figure [3j 

In section [H we consider the behaviour of our solution of the equation of motion obtained 
in section [2] in the long-time limit. The solution appears to be incompatible with a statistically 
stationary limit, but this is shown not to be the case. We also show that the probability distribution 
of the gradient of the angle has an approximately log-normal distribution. This is consistent with 
the existence of apparent singularities in the rod textures, where the angle of the rods changes very 
abruptly. Section [5] summarises the results and discusses how the patterns observed at long times 
can be understood. The numerical simulations are described in an appendix. 
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Fig. 4. — Given a non-oriented vector field n(r) in two dimensions and a closed curve C, the 
Poincare index N{C) is defined as the number of multiples of 27r by which the direction of n rotates 
(in the clockwise direction) as C is traversed (also clockwise). For a non-oriented vector field, 
such as the direction of the rods, the Poincare index may take half-integer values, (a) For a field 
without singularities, n = 0. (b) For a curve which encircles a core, iV = i. This singularity can 
be regarded as having a charge of N = ^. (c) For a curve which encircles a delta, = — ^. (d) For 
a curve which encircles more than one singularity, their charges are summed. This curve encircles 
a combination of a core and delta which is termed a loop. For this case A^ = — ^ + ^ = 0. 



2. Equation of motion and its solution 
2.1. Derivation of the equation of motion 



While the equation of motion which we consider is a limiting case of that derived by iJefferv 



(|1922| ). the general calculation is quite lengthy and insight is gained from a simple derivation. 
Strictly speaking, in the calculation below we consider the motion of dumbells, that is pairs of 
particles (which are dragged by the fluid) connected by a rigid rod (which is not influenced by the 
flow). However, the equation of motion we obtain is independent of the length a of the rod in the 
limit as a — > 0, and by imagining a rod as being formed by overlaying dumbells of different lengths, 
we surmise that our equation describes a short symmetric rod with a general distribution of its 
viscous drag along its length. The rods are advected by a velocity field v{r,t), which is characterised 
by a correlation time r, correlation length ^ and typical magnitude vq. In a multiscale turbulent 
flow, it is the correlation time and correlation length of the smallest eddies which are relevant 
here (that is, we identify r and with the Kolmogorov time and the Kolmogorov length of the 
turbulence, respectively). 



Fig. 5. — The direction vector n (black lines) is asymptotic to the vector field of eigenvector field 
11+ (red lines) . The vector field it_|_ is undefined in gyres, where the normal form of the monodromy 
matrix is a rotation. The Poincare index of the field 1x4. on the boundary of the gyre need not be 
equal to zero, whereas the Poicare index of n is zero. In these cases the field n rotates by tt in the 
vicinity of one or more scar lines (green). 

The configuration of the rod can be specified by the position r{t) of its mid-point at time t, and 
by a unit vector n{t) aligned with the rod (the binary ambiguity of the evolution this is resolved 
by requiring continuity). The rods have an initial direction ng, which is a smooth function of the 
position r. Our aim is to obtain equations of motion for r and n, using these to understand the 
vector field n{r,t) describing the orientation of the rods which have reached position r at time t. 

We simplify by assuming that the rod length a is very short compared to the correlation length 
^: a/(^ <^ 1. The centre of the rod with position r is therefore assumed to move according to the 
advective equation of motion, r = v(r{t),t) (we neglect small O^a^) corrections). To obtain the 
equation of motion for the direction of the rod, we use a linear approximation for the velocity 
difference 5v between the centre of the rod, r, and one of the particles at its ends, at r + 6r: 



where A{r,t) is the strain-rate matrix (a 2 x 2 matrix with elements Aij = dvi/drj, which satisfies 
tr[A] = because V • = 0). The line between the two particles has direction specified by the unit 
vector n. The equation for force balance on one of the particles at the end of the rod is r = v — Tn, 
where v and r are evaluated at the position of the particle at the end of the rod and where T is 
proportional to the tension in the rod, which keeps the separation of the two particles at its ends 
constant. Because r = d at the centre of the rod, we obtain 6r = 6v — Tn, and the equation of the 
constraint is 6r ■ n = 0. From these we find T = 6v ■ n. Combining these results with ([2]) we find 
an equation of motion for n: 



6v{r, t) = A(r, t)6r 



(2) 



h = An — (n • An)n . 



(3) 
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Fig. 6. — (a) The rod direction (black) is a smooth vector field containing a scar line which ends 
on the boundary of a gyre. As time increases, the scar line narrows (b). When the scar line has 
narrowed to the extent that it does not include the actual position of any rod, it disappears (c). 
This leaves a point singularity at the end of the scar line: in this case a delta. In practice, the 
picture is more complex because the positions of the gyre and the scar line both change as time 
increases. 



This equation of motion is the same as that obtained by iJefFervl (|l922l ) for a prolate ellipsoid of 
rotation, in the limit as the aspect ratio approaches infinity. 



2.2. Solution of the equation of motion 

We now consider how a solution of the equation of motion (I3D may be obtained from the mon- 
odromy matrix of the flow. The monodromy matrix M describes the evolution of the infinitesimal 
separation vector br of two points advected by the flow, r = v{r,t): we write the separation of two 
points at time t in the form 

6r{t)=M{r{t),t,to)6r{to) . (4) 

Note that M is written as a function of the position r reached by the rod at time t, and of the final 
and initial times, t, to, respectively. The monodromy matrix satisfies the differential equation 

^M = A{r{t),t)M (5) 
at 

where r{t) is the trajectory of the centre of the rod. The initial condition for equation ([5]) is 
M(r,toiio) = I; where I is the identity matrix, for all positions r. Now define no(r'o) as the initial 
direction, at time to, of the rod, expressed as a function of the initial position Tq. Let us consider 
the vector field 

a{t) = M{r,t,to)no{ro) (6) 

where rQ{r,t,to) is the initial position, at time tQ, of a rod which reaches r at time t. If we write 
a(t) = a{t)n{t), we find that n(t) satisfies the equation of motion ([3]) above. Also, it satisfies 
the initial conditions: n(to) = no(r(to), *o)) since M.{r,to,tQ) = I. Thus we can determine the 
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Fig. 7. — Illustrating the correspondence between the rod direction field n(r,t) (black) and the 
eigenvector field u^{r,t) (red) at large time. 



orientation vector of the rods from the monodromy matrix by normalising the vector a{t): 

M(r,t,to)no(ro) 
^ |M(r,t,to)no(ro)| 

(where the initial position ro is a function of r, t, Iq). 



(7) 



2.3. Asymptotic form of the solution 

Let A+ and u+ be respectively the eigenvalue of M with the largest magnitude and the cor- 
responding eigenvector, normalised to unit length. We term these the dominant eigenvalue and 
eigenvector. The other eigenvalue and eigenvector are denoted A_ and u and are termed subdom- 
inant. We assume that the random flow has a positive Lyapunov exponent, so that the elements 
of M(r,t, io) tend to increase exponentially as a function of t — to- Correspondingly |A+/A_| is 
expected to increase exponentially (recall that A+A_ = 1). If |A+/A_| » 1, applying the matrix M 
to almost any vector is expected to result in a vector which is nearly aligned with u+. In particular, 
as t — to ^ CO we expect that n(r, t) ~ u^{r,t) for almost all points in the plane. This is illustrated 
by the simulation in figure [71 
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3. Apparent singularities of the direction field 

Here we consider whether it is possible for the vector field n(r,t) to have singularities, where 
n changes discontinuously as a function of r. First we show (section 13. ip that it is not possible for 
n{r,t) to have singularities in a strict sense. It is however possible that the field could approach 
a singularity in some asymptotic sense. Accordingly, we also consider (section I3.2p whether the 
eigenvector field u^{r,t), to which n(r,t) is asymptotic, has any singularities. Although u^{r,t) 
does not have singularities, we show that it can have a non-trivial topology. There are regions 
where the monodromy matrix M is elliptic (with conjugate eigenvalues on the unit circle) so that 
the dominant eigenvector it^- is not defined. We term these regions of rotational flow gyres. We 
find that the Poincare index of the eigenvector it+ around the boundary of a gyre can be non-zero. 
In sections 13.31 13.41 we consider how the smooth field n{r,t) can be asymptotic to the topologically 
non-trivial field u^{r,t). 

3.1. Absence of singularities 

The monodromy matrix M(r, t, to) is a smooth function of the final position of the trajectory, 
r. The solution ([7]) can therefore only be discontinuous if the initial direction field is discontinuous, 
or if the denominator |Mno| is equal to zero, which is only possible if there are points where 
det(M) = 0. This is not possible since we consider area-preserving fiows, where det(M) = 1. If 
the initial direction vector field no(r) is non-singular, we therefore conclude that the direction field 
n(r,t) remains non-singular for all times. Because the vector field generated by ([7|) is smooth, the 
Poincare index of this field is zero for any closed curve, in apparent contradiction to the simulations 
shown in figure [TJ 

3.2. Topology of the eigenvector field 

We have shown that the direction field n(r, t) is asymptotic to the field of eigenvectors, u+{r, t). 
We shall see that the latter field has a non-trivial topology. 

The only type of singularity of the eigenvector field which is possible is where the monodromy 
matrix is equal to the identity matrix. It is a co-dimension three condition for the monodromy 
matrix to have this form, so it is non-generic in the two-dimensional problem which we consider. 
There is, however, another way in which the eigenvector field u^{r, t) can have non-trivial topology. 

In an area-preserving fiow there will be regions of the plane where the eigenvalues are complex 
and have the same magnitude. We refer to these regions where the normal form of M(r,t) is a 
rotation as gyres. Each gyre is surrounded by a boundary. We find that the Poincare index of the 
field w+ on the boundary of a gyre may not be equal to zero (two examples are illustrated in figure 
[8|). This appears to contradict the result that n is asymptotic to ti+, because we have seen that 



Fig. 8. — Eigenvector field w+(r, t) (red) and rod direction field n(r, t) (black) at small time t — to- 
The regions where there are no red vectors arise because the eigenvector field of the dominant 
eigenvalue of the monodromy matrix M(r,t,to) is undefined in regions (which we term gyres) 
where M(r,t,to) is elliptic. Note that at small t — to the gyres occupy a large fraction of the area 
and have a simple boundaries. In this figure it can be seen that two of the gyres have non-zero 
Poincare index. 

the Poincare index of n is always zero. 

3.3. Asymptotic singularities of the direction field 

We have seen that n{r,t) is non-singular, but that it is asymptotic to a vector field u^{r,t) 
which may be topologically non-trivial. One way to resolve this contradiction is to assume that 
the field u+(r,t) has become trivial by the time n(r,t) approaches it, due to gyres with opposite 
topological charges coalescing. 

There is, however, another route to resolving this apparent contradiction which is both more 
interesting and which does lead to an explanation of the textures seen in figure [TJ Let us consider 
the set of points where n need not be asymptotic to u+. We write the initial direction field as 



no = a+w+ -|- q;_u_ 



(8) 
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The vector n is proportional to Mno = a+A+w^ + a_A_M_. The ratio of eigenvalues, A+A- 
exp(27|t — tol) grows exponentially, with Lyapunov exponent 7, as t — to increases. Hence n aligns 
increasingly closely with except when is sufficiently small. The locus where a+ = forms a 
set of lines in the plane, and as we cross these lines the direction of n rotates by itTr. We term these 
lines scar lines. The vector n differs significantly from when |a_|_| exp(27|t — tol) = 0(1). This 
region where the direction flips therefore becomes vanishingly small at t — to ^ 00. Accordingly, 
we can think of the scars lines as 'healing over', that is, becoming invisible. 

The scar lines must terminate at gyres. Figure [5] is a schematic illustration the fields n and 
in the vicinity of a charged gyre and its associated scar line. Figure [9] shows scar line in our 
numerical simulations. 

We conclude this section by remarking that the eigenvectors w+ and tt_ become co-linear 
on the boundary of the gyre. This observation can be understood using the following argument. 
On the boundary of the gyre, the matrix M only has one eigenvalue (which may be +1 or — 1). 
The set of 2 x 2 matrices satisfying det M = 1 has three parameters, and if the eigenvalues are 
constrained to be A = 1 (say), it becomes a two-parameter family of matrices. We now identify 
a parametrisation of this two parameter family. Consider the eigenvalue equation, Fit = Xu, for 
matrices of the Jordan form 

These are a one parameter family of matrices which have only one eigenvector, u = (1,0), and one 
eigenvalue, A = 1. If R is a rotation matrix 

f cosO si„9\ 
\ — sm d cos t) J 

we see that we can generate a two-parameter family of 2 x 2 matrices M(^, k) = R^^(0)F(k;)R(0) 
which have only one eigenvalue, A = 1. By construction of the matrix M(0,k) we can show that 
this two parameter family spans the set of 2 x 2 matrices with only one eigenvalue, A = 1. But 
we have seen that these matrices have only one eigenvector, namely u = R(^)(l, 0)^. We conclude 
that as we approach the boundary of a gyre from the outside, the two eigenvectors and U- 
become co-linear. This implies that a+ and a_ both diverge as we approach the boundary of the 
gyre. 

3.4. Disappearance of scar lines and emergence of point singularities 

As noted in section [3T3l above, the width of the region around a scar line where the fields n and 
M-i- are significantly misaligned shrinks as t — to ^ 00. As this region shrinks, eventually there is 
a small probability that any rod actually lies in the region where these vectors are misaligned. In 
this case, for all practical purposes the scar line has disappeared. Consider a loop which encircles 
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Fig. 9. — Numerical examples of scar lines. The rods directions n are shown in black, eigenvector 
ti+ is shown in red, and the position of the scar line is indicated by a sampling of points where 
< 10^^ (green). The direction of the rods is seen to flip by around in the vicinity of the scar 

line. 

the end of a scar line. Initially the Poincare index of n about this loop is zero. When the angle 
change of zbvr associated with crossing the scar line disappears, the Poincare index of the circuit 
becomes = it^. The disappearance of the scar line is therefore associated with the emergence 
of a point singularity at the positions where the ends of this line were located. This is illustrated 
schematically in figure [H and by the numerical simulations in figure [TOj This effect gives rise to 
the apparent singularities seen in figured! 



4. The long-time limit 

4.1. Sensitivity to final position 

If we assume that the velocity field v{r,t) is statistically stationary, we expect that at long 
time the patterns formed by the rods also become statistically stationary (so that at long times it 
becomes impossible to estimate the time from the statistics of a realisation of the rod positions). 
This property is, however, not manifest in the solution d?]). As t — Iq ^ oo, the norm of the 
monodromy matrix grows. Also, while its elements are everywhere a smooth function of the initial 
position, the elements of M(r, t, to) do become ever more sensitive to the position r as t — to ^ oo. 
These observations suggest that as time increases the vector field n{r,t) should vary increasingly 
rapidly as a function of r, the final position of the rods. We argue below that this is not the case, 
and that n(r,t) does reach a statistically stationary state. 

If the eigenvector u^{r,t) corresponding to the largest eigenvalue is less sensitive to the final 
position r than the matrix M(r,t, to) itself, then the apparent contradiction discussed above can 
be resolved. We now argue that this is, in fact, the case. 
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Fig. 10. — Simulation showing healed scar lines. The rods are shown in black, and the eigenvector 
field is shown in red. The points in green show positions of rods where |a+| < 10~^ (and are 
therefore very close to scar lines). In these examples, because the scar line has become very 
narrow, the direction of the rods is not seen to flip around in the vicinity of the scar line, and the 
end of the scar line is marked by an apparent point singularity, namely a loop (a) or a delta (b). 

Let us consider a 2 x 2 random matrix M(i, to) generated by an equation of motion M = A(t)M 
(that is, by equation ([5])), where A(t) is a traceless 2x2 matrix generated by a stationary random 
process. We apply the initial condition M(to,to) = I; where I is the unit matrix. The time- 
dependence of the matrix M has a positive Lyapunov exponent, 7, describing exponential growth of 
the largest eigenvalue A+ as a function of |t — to|- Our discussion of the sensitivity of the solution will 
use an observation about the subdominat eigenvector w_, corresponding to the smallest eigenvalue 
A_. We start by showing that this eigenvector approaches a constant direction, which depends upon 
the initial realisation of A in the first few multiples of the Lyapunov time, 7~^. (The direction of 
the other eigenvector, continues to fluctuate randomly.) 

To demonstrate this result, we consider the change in this eigenvector U-(t) during a small 
timestep 6t. Let X(t) be a matrix formed using the eigenvectors of M(t), such that D = XMX~^ = 
diag(A+, A_) is the diagonal matrix formed from the eigenvalues of M(t). In the transformed basis, 
the matrix at time t + 6t is 

M' = XM(t + 6t)X-^ = X[I + A(^t]MX-i = [I + A'{t)6t]B (11) 

where A' = XAX^^. We write the eigenvalue equation for the subdominant eigenvalue of the 
matrix M, with subdominant eigenvector u^: 



M'u'_ 



1 + A[^6t A'^^St 
A'^iSt 1 + A'226t 

X+{l + A[^6t) X^A[26t 
X+A'^.S A_(l + A'225t) 





(A_ +a_) - . (12 
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Neglecting terms of higher order in 6t, the first element of this eigenvalue equation gives 

5u'_ = -^-^-^A[^6t . (13) 

In the limit as t ^ oo the eigenvalues satisfy |A+/A_| oo, and 5u'_/5t 0. We therefore conclude 
that the eigenvector of the subdominant eigenvector approaches a constant direction. Writing the 
eigenvector of the dominant eigenvalue of M' as it+ = (1, Su'_^_), the corresponding expression is 

Su+ = A'2i6t . (14) 

Here the coefficient of 5t does not approach zero as t — > (X), and we conclude that the dominant 
eigenvector continues to rotate in the large time limit. 

Now given the orientations of the rods at time t, consider their orientations at the earlier time 
to- This map is determined by a time-reversed version of equation Its solution is constructed by 
analogy with ([7]), replacing M with M~^. The eigenvector of M(t,to) corresponding to its largest 
eigenvalue is also the eigenvector of M^^(t,to) corresponding to its smallest eigenvalue. Using 
the result discussed above, the eigenvector corresponding to the smallest eigenvalue of M^^(t,to) 
becomes insensitive to to when — tol ^ 1- Correspondingly, the eigenvector of M(t,to) corre- 
sponding to the largest eigenvalue becomes insensitive to to- We conclude that although the matrix 
M(t, to) has an increasingly sensitive dependence upon position as t — to — >■ oo, the eigenvector ti+ 
does not become increasingly sensitive. Because the rod directions are asymptotic to these vectors, 
the rod directions do not become increasingly sensitive to the position r as time increases. 

There are regions where the matrix M is not hyperbolic, so that there is no largest eigenvalue 
and consequently it+ is not defined. However, as |t — to| oo, the norm of M increases almost 
everywhere, and the fraction of the area of the plane occupied by regions where ti+ is not defined 
approaches zero. 

We conclude that as t — to ^ oo, the vector field n(r, t) is statistically stationary, approaching 
the vector field u+(r,t). 



4.2. Distribution of angle gradients 

We have seen that the rod directions do not become increasingly sensitive to position as 
time increases. It is desirable to quantify the sensitivity to position. We have seen that the 
textures formed by the rod orientations show regions where the rod direction varies very rapidly 
with position, relative to other regions. Earlier, we described how the existence of scar lines 
explains the structures seen in specific realisations of the patterns. In this section we consider the 
probability distribution of the angle gradient, showing that the distribution is very broad, being 
well approximated by a log-normal distribution. This very broad distibution of the angle gradient 
is consistent with the existence of the structures described in section [3l 
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We now consider how to calculate the angle gradient g = V0. In the following, we obtain 
an expression for one component, gi, g. We obtain an equation for gi, equation (jlSp . which 
is easily argued to be log-normally distributed. It is, however, less clear that this formula for gi 
gives results which are well defined. We discuss this point in some detail after deriving (jlSp . before 
finally presenting a brief argument that gi is approximately log-normally distributed at the end of 
this section. 

Consider the difference between the eigenvector direction between two monodromy matrices 
evaluated along neighbouring trajectories. The reference trajectory has monodromy matrix M(t) 
and the neighbouring trajectory has monodromy matrix M(t) + (5M(t). We have seen that the 
subdominant eigenvector it_ of each monodromy matrix approaches a constant direction as t — > 
oo, so the angle between them, S9{t), must approach a constant value, that is 69(t) 69oo as 
t oo. Let 5M(t) be the change in the monodromy matrix due to shifting the end point of the 
trajectory at time t from r = (ri,r2) to r + 6r = (ri,r2) + (5ri,0). The first component of g is 
gi = limsr^^oSO/Sri. 

We introduce an orthonormal basis ui, U2 satisfying Wj • uj = 5ij, where U2 = U-{t) is the 
subdominant eigenvector of M(i). The elements of M in this basis are M-j = Ui ■ 'M.{t)uj, which 
form the matrix 

When the end-point of the rod trajectory is shifted by a distance 6r = {6ri,0), the matrix M' 
is perturbed to M' + 5M.' , and the angle of the subdominant eigenvector ii_ changes by a small 
amount 69, which can be obtained by solving the eigenvalue equation 



'22 




ML+5ML X_ + 6ML I \ 1 I ' M 1 



(A_+<5A_) ^ . (16) 



Using the first line of this equation to solve for 69, retaining leading order terms we obtain 

- - (17) 



- A 



Note that when t is large, so that A+/A_ ^ 1, we may drop the term A_ from the denominator, 
and approximate the first element of the gradient vector by 

dM(t) 

gi = lim — ~ ^ ^, , . (18) 

sr^^Q Sn Ul ■ m{t)ui ^ ' 

The angle gradient must approach a definite value as t — > 00, but it is not immediately clear that 
this expression approaches a constant value. We must look at (jlSp more carefully to see why this 
is in fact true. 
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It is desirable to have an explicit expression for the coefficients 5M[-{t). Note that the mon- 
odromy matrix M(t) satisfying dM/dt = A(f)M can be approximated by a product: 

Int(t/<5t) 

M(t)= hm n [l + Am] . (19) 

i=i 

Writing B = dA/dri, the monodromy matrix for the displaced trajectory is 

Int(t/5t) 

M{t) + SM{t) = lim Yl [1 + A(j5t)5t + B(j5t)5ri5t] 

i=i 

lnt{t/5t) lnt{{t-t')/St) 

= 5ri lim Yj n + + •^''^*)'^*] 

k=i j=i 

Int {t' /St) 

xB{k6t)6t Yl [^ + Mj^t)]+0{B^) . (20) 



We find 



dM{t) 
dri 



[ dt' M{t,t')B{t')M{t',0) . (21) 



Having obtained an expression for 5M, we return to considering why gi, given by equation 
(jlSp . is independent of t in the limit as t — s- cxd. Let us introduce the initial time in the arguments 
of the monodromy matrix, writing the monodromy matrix giving displacements at time t in terms 
of those at time to as M(t,to)- Consider the vectors vi = M(t, 0)iti and V2 = (5M(t, 0)w2, where 
til, U2 are two arbitrary vectors. We will show that the vectors i^i V2 almost always become 
co-linear as t — > oo. First choose a time ti such that (t — ti)^ ^ 1. Note that we can write 
M.{t,to) = M(t, ti)M(ti, to)- The direction of the vector vi = M.{t,to)ui is almost always nearly 
co-linear with the direction of the dominant eigenvector of M(t, ti), independent of the vector ui. 
On the case of the vector V2 = S'M.{t,tQ)u2, note that we can write 



6M{t,to) = 6riM{t,ti) [ ' dt' M{ti,t')B{t')M{t',to) 

+Sri [ dt' M{t,t')B{t')M{t',to) 
Jti 

= M{t,ti)6M{t,,to)[l + 0{{t-ti)/t)] (22) 

so that to leading order V2 is also co-linear with the dominant eigenvector of M(t, ti). We conclude 
that the vectors i^i and V2 are almost always co-linear, provided ^{t — to) ^ 1. 

Let us consider the evaluation of (llSp in the case where 



M = M(t,to) = M(t,t')M(t',to) = M2M1 



(23) 
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where Mi = M(t',to)) M2 = M(t,t',). Correspondingly, neglecting terms of order 6rf, we have 

6M = M2(5Mi + ^MaMi . (24) 

We consider the case where 7(^ — ^0) ^ with t > t' > tQ. In order to establish that the angle 
69 becomes asymptotically independent of time, we must show that 69 = 69i, where 69 is given by 
()18|) and where 69i is the expression obtained by replacing (^M, M with 5Mi, Mi. Thus (in view 
of ([23|) and (p^ ) we must demonstrate that 

ui-6MiU2 ui-M26MiU2 ui ■ 6M2M1U2 

69 = = . (25) 

111 • Mi^i Ml • M2M1U1 ui ■ M2Miiii ^ ^ 

The second term on the right-hand side of the equality is negligible, because M1U2 = A_u_, and 
A_ — > as t — > 00. In the first term the additional factor of M2 makes no difference to the value 
69 only if the vectors 6M.1U2 and MiUi are co- linear. But we argued above that these vectors are 
asymptotically co-linear in the limit as 7(t — to) ~^ 00. Thus we conclude that the angle 69 between 
two sub-dominant eigenvectors ii_ in forward-time propagation does become independent of time 
as t — > 00, justifying p^ . We can now use the arguments of section [d] to draw conclusions about 
the dependence of the reverse-time propagation of the dominant eigenvectors u+, which determine 
the rod direction. We conclude that the angle gradient at time t does become independent of the 
initial time to as t — to — > 00. 

We have seen that the angle gradient of the orientation field of the rods remains finite in 
the long-time limit. It is of interest to consider the probability distribution of the angle gradient. 
We shall argue that this quantity has an approximately log-normal distribution. We note that 
the monodromy matrix may be expressed as a product of a large number of independent random 
factors. It is clear that when 7t ^ 1, the distribution of the matrix elements of both 6M. and M are 
log-normal. The distribution of their ratio is also log-normal. We conclude that the distribution of 
the angle gradient, given by (jlSh . is therefore also log-normal at large times, provided the correlation 
time T is short compared to the Lyapunov time 7~^. In figure [TT] this result is illustrated by a 
histogram of the distribution of the logarithm of the angle gradient for the same parameter values 
as used in the other numerical simulations. A Gaussian fit matches the histogram very closely. 



5. Discussion of the rod textures in the long time limit 

We have shown that the rod textures seen on figured] may be understood in terms of concepts 
introduced in section [3l We showed that the direction field n(r,t) is asymptotic to the vector field 
of the dominant eigenvector u^{r, t) of the monodromy matrix M(r, t). However, we observed that 
n has a simple topology, whereas w+ has a non-zero Poincare index upon traversing boundaries 
of some of its gyres. In order to reconcile the different topologies of these fields, we noted that 
this asymptotic correspondence breaks down on scar lines, where the direction vector n rotates 
abruptly by vr. We showed that the size of the region where the direction reverses can decrease 
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Fig. 11. — Histogram of the probability density of the angle gradient, showing that this has an 
approximately log-normal distribution. The parameter values of the simulation are the same as for 
the other figures. 

as the norm of the monodromy matrix increases, so that these scar lines can heal over when they 
become sufficiently narrow that it is unlikely that a rod lies in the region of the scar line. When 
the scar line has healed, there appears to be a point singularity with non-trivial topology at each 
of its ends. 

In the long-time limit, the application of these concepts becomes increasingly problematic. This 
is because, at very large times, the monodromy matrix M(r, t) becomes increasingly sensitive to the 
final position of the rods, r. As t — to increases the gyres may shrink in area, their boundaries may 
stretch, and they may merge together. Also, gyres with a non-zero Poincare index may disappear. 
We expect that at very large times the gyres are extended into lines where tr(M) changes sign. 
These lines are expected to become ever more closely spaced as t — to ^ oo, with typical spacing 
^exp(— A|t — to\). Following the reasoning presented in section [3l at large t — Iq we expect that 

fluctuates on a lengthscale ^, independent of t — to- In this limit most of the gyres must have 
Poincare index equal to zero. 

The trajectory of the scar lines, which depends on both eigenvectors, also becomes an increas- 
ingly sensitive function of r as t — to increases, until the scar line is densely folded but its region 
of effect is exponentially narrow. As well as healing over, the scar lines may also stretch and fold, 
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resulting in rods with an orientation which differs from that of the surrounding rods, apparently 
randomly scattered in the plane. Examples of these erratically misaligned rods can be seen in figure 

m 

The reason for these pathologies is that the definition of the gyres and scar lines refers to the 
initial conditions of the problem (specifically, the initial direction field no), which become irrelevant 
in the long time limit. 

What can we say about the rod textures at long times? We have shown that the patterns are 
statistically stationary for long times, because the direction of the dominant eigenvector M_|_(r,i) 
is determined only by the recent history of the monodromy matrix, over a few multiples of its 
Lyapunov time, 7~^. We have also seen that the distribution of angle gradients is approximately 
log-normal. This very broad distribution is consistent with the apparent singularities which we have 
discussed in section [3l The patterns which are seen at a very large time t can be understood by 
applying the same principles as are used to understand the patterns at short times. We assume that 
the direction field n(r, to) at time is known, where t — to is of order the Lyapunov time, 7~^. We 
know that the direction field no(r, to) is smooth, although it may have apparent singularities of the 
types discussed in section [3j In the time between to and t, this field will undergo further evolution 
involving the production of additional apparent singularities, which can be analysed by considering 
the gyres associated with the monodromy matrix M(r, t, to) and the scar lines associated with the 
initial orientation field no(r,to). At the same time, the apparent singularities which are already 
present in the initial orientation field n(r,to) become less visible as the pattern is stretched and 
folded. We conclude that the analysis of section [3] is sufficient to explain the nature of the textures 
seen at large times. 
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7. Appendix 



Numerical simulations used a synthetic vector field v{x, y, t) which was periodic in x, y (with 
period L) and in t (with period T), generated from a random stream function il^{x,y,t): the 
components of the velocity field are Vx = dif^/dy, Vy = —dijj/dx. The stream function is written in 
terms of its Fourier decomposition 



,Pix,y,t) = 5;^^A(fc„fc„a;)e<*^--+'=-^+-*) 



(26) 



where kx, ky are integer multiples of 2n/L and where u is an integer multiple of 2'k/T. The 
coefficients A{kx, ky,u) are random Gaussian variables with the following properties 

{A{kx,ky,Lu)) = 
{A{kx,ky,u;)A*ik',,k'^,u;')) = 6k^k'Jkyk'J..'ivoO\2^?^'j^e^P + 



A(^kx,ky,Lo^ — A { kx, ky, lo^ . 
The correlation function of il){x,y,t) is given by 

{x - x'f + {y- y'f 



{il:{x,y,t)ij{x',y',t')) = (i^oO exp 



2e 



exp 



l\2 



it - t') 

2r2 



(27) 



(28) 



This preprint was prepared with the AAS IM^^X macros v5.2. 
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The fast Fourier transform was used to calculate Fourier components at discrete time steps t„ = n5t. 
In the simulations we used r = 0.1, ^ = 0.1, vq = 1.0. 

In all of the simulations the rods were all initially in the same direction, that is uq was 
independent of r. 

The colour mapping of figure [2] was produced using MATLAB. 



